Will Nonlinear Peculiar Velocity and Inhomogeneous Reionization Spoil 21cm 
Cosmology from the Epoch of Reionization? 
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The 21cm background from the epoch of reionization is a promising cosmological probe: Hne-of- 
sight velocity fluctuations distort redshift, so brightness fluctuations in Fourier space depend upon 
angle, which linear theory shows can separate cosmological from astrophysical information. Nonhn- 
ear fluctuations in ionization, density and velocity change this, however. The validity and accuracy 
of the separation scheme are tested here for the first time, by detailed reionization simulations. The 
scheme works reasonably well early in reionization (< 40% ionized), but not late (> 80% ionized). 

PACS numbers: 98.80.Bp,98.58.Ge,95.75.Pq 



Introduction. Neutral hydrogen atoms in the inter- 
galactic medium (IGM) at high redshift contribute a dif- 
fuse background of redshifted 21cm-line radiation which 
encodes a wealth of information about physical condi- 
tions in the early universe at z > 6, during and before 
the epoch of reionization (EOR). To derive cosmological 
information from this, however, we must be able to sep- 
arate the dependence of the signal on the cosmological 
properties of the background universe (i.e. total matter 
density fluctuations) from that on the complex astrophys- 
ical processes that cause the thermal and ionization state 
of the intergalactic gas to fluctuate (e.g. energy-release 
from star and galaxy formation). The anisotropy intro- 
duced by the peculiar velocity of the gas, induced by 
structure formation, is the key to this separation. 

According to linear perturbation theory, the power 
spectrum of the 21cm background fluctuations can be 
expressed as a sum of terms which depend on different 
powers of the cosine, /ik, of the angle between the line- 
of-sight (LOS) n and the wavevector of a given Fourier 
mode k [1]. Different terms represent contributions from 
different sources of fluctuations, including fluctuations in 
the total matter density, velocity, and hydrogen ionized 
fraction, and thereby in principle provide a means of sep- 
arating the effects of cosmology and astrophysics. In par- 
ticular, future measurements, it is proposed [l|, can be 
used to fit this theoretical dependence of the power spec- 
trum on /ik to extract the power spectrum of total matter 
density fluctuations ~ the cosmological jewel. 

The success of this approach, however, depends upon 
the validity of the linear /ik-decomposition. While it 
is true that fluctuations in the total matter density at 
high redshift are likely to be of linear amplitude on 



large scales, the nonlinearity of small-scale structure in 
density, velocity and reionization patchiness can leave 
its imprint on the 21cm signal, which might result in 
nonlinear distortion in the 3D power spectrum of 21cm 
brightness temperature fluctuations and so spoil the lin- 
ear /ik-decomposition. In what follows, we will exam- 
ine this question. We will assess the accuracy of this 
method for deriving cosmological information from the 
21cm background by using the results of new large-scale 
N-body-|-radiative transfer simulations of cosmic reion- 
ization as mock data from which to "measure" the mat- 
ter density power spectrum. Our simulation volume 
(425 Mpc//i)'^ is large enough to make the sampling errors 
smaller than the systematic errors, while approaching the 
size of upcoming surveys like LOFArHi (5° X 5°). 

The 21 cm brightness temperature in redshift 
space. The observed frequency fobs reflects both 
the cosmological redshift Zcos from the time of emission 
and the Doppler shift associated with the peculiar ra- 
dial velocity there, i.e. fobs = fo/(l + ^obs), where 
1 + Zohs = (1 + Zcos)(l ~ ~^)~^- In the "distorted" co- 
moving coordinate system known as redshift space, the 
position of the emitter is the apparent comoving position 
if the redshift is interpreted as cosmological only, which 
shifts the real comoving coordinate r along the LOS to 
s = r + (f II /ai7(2)) n, where a = (1 + Zcos) "'^ is the 
cosmic scale factor. Henceforth, superscripts "r" and "s" 
denote quantities in real- and redshift-space, respectively, 
and we will write z for Zcos unless otherwise noted. 

Peculiar velocity can enhance/suppress the 21cm 
brightness temperature measured relative to that of the 
cosmic microwave background (CMB), from over/under- 
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Here we focus on the hmit where the spin temper- 
ature TJ ^ TcMB, vahd soon after reionization be- 
gins. As such, we can neglect the dependence on spin 
temperature, but our discussion can be readily gen- 
eralized to the case of finite TJ. In equation ([T|), 

ST.iz) = (23.88mK)(t^)^ 



ih^ 10 -^HI,: 



{z) is the 



mean 21cm signal (where xiii,m{z) is the mass-weighted 
global neutral fraction), Sp^^{r) is the neutral hydrogen 

overdensity, and Sg^^{r) = -g^^(r) is the gradient of 
proper radial peculiar velocity along the LOS, normalized 
by where H{z) is the Hubble parameter. 

The 21 cm power spectrum anisotropy in redshift- 
space. In Fourier space, the three-dimensional 
power spectrum of 21cm brightness temperature 
(hereafter "21cm power spectrum") is defined as 

(Wt* {k)Wb{k')) = (27r)3F^j,(k)5(3)(]^_k'). There are 
three approaches to model the 21cm signal in redshift 
space, based on different assumptions on peculiar veloc- 
ity and reionization patchiness, as follows. 

• Linear scheme (linear velocity- density relation, lin- 
earized neutral fraction fluctuations). As described 
above, the linear scheme was originally proposed in the 
context of linear perturbation theory [l| and later re- 
derived with weaker assumptions Q which are: (1) the 
velocity (v'') and total density fluctuations (Jp satisfy 

Hiz)' 



the linear relation, z;f|'(k) = i (^i+jj ^p(k)/^k/fc, (2) the 
baryon distribution traces the CDM, (3) the peculiar ve- 
locity, the hydrogen density fluctuation (Sp^^), and the 
neutral fraction fluctuation {S^^^ ) are all linearized. Un- 
der these assumptions, the 3D 21cm power spectrum can 
be expanded in polynomials of /ik = k • n/|k|, 

P^'^"(k,z) = F^o(fc,z) + P^2(fc,z)A*^ + P^4(fc,z)4, (2) 

where the moments {P^^o , Pp2 , P^4 } are functions of only 
k = |k| and redshift z. More importantly, the 4**^- 
moment depends only on the total density power spec- 
trum pj'*°*^' and the global neutral fraction XHi,in('2^), 



p^.(fc,z)^<5r-(z)p;'*°'^'(fc,z) 



(3) 



while other moments are "contaminated" by power spec- 
tra due to reionization and/or spin temperature. In prin- 
ciple, then, cosmological information can be extracted 
from the 21cm signal by fitting the measured P^rp(\i,z) 
to equation 1^ to isolate the 4'^ moment. 

• Quasi-linear fi]i^- decomposition scheme (linear velocity- 
density relation, linearized neutral overdensity). The 



assumption of 5^^^ ^ 1 breaks down when a;Hi,m > 0.5, 
which results in a significant contribution to the power 
spectrum from higher-order terms, e.g. PJ i a lii, 
etc. Fortunately, the 3D 21cm power spectrum can still 
be written in the form of equation ([2]), and, more im- 
portantly, the 4*'^-moment is still related to the matter 
power spectrum as in equation ([3]), if we adopt the same 
assumptions (1) and (2) as in the linear scheme, but as- 
sume peculiar velocity and the neutral density fluctua- 



tion SI, 



^^Hi + ^ph'^^hi' as opposed to S^^ 



alone, are linearized. This is the so-called quasi-linear 
/ik-decomposition formalism Q , in which only lower mo- 
ments differ from the linear scheme prediction. 

• Fully-nonlinear scheme (nonlinear velocity, nonlinear 
neutral overdensity). In the optically-thin approxima- 
tion, two nonlinear effects of peculiar velocity must be 
taken into account: (1) the 21cm brightness temperature 
is corrected for velocity gradient as in equation ([1]); (2) 
when the rea/-space comoving coordinates r are mapped 
to redshift-space coordinates, volume elements are also 
resized according to 6V^{s) — SV^{r) |l -I- (5g^^(r)|. For- 
tunately, the combined effect allows us to compute the 
21cm brightness temperature in redshift-space with a 
simple formula [3] 



(4) 



where (5^^^ (s) = ?1hi(s) /niii{z) — 1 is the neutral overden- 
sity in redshift-space. 

Angular separation of 3D power spectrum. Given 
the observed 3D power spectrum on a uniform 3D grid in 
k-space, its moments can be decomposed using the x^-fit 
as follows. For a given LOS, 3D modes P^j,(k) with the 
same fi^ and same k but different azimuthal angle are av- 
eraged to give a measure of P^rp{k, /ik) and the associated 

sampling variance ap{k,fi]^) = Plr^l^'Mk), 
where N^^^^k is the number of modes with the same /ik 
and k. For the 3D grid from the simulated data de- 
scribed below, we further combine the measures along 
three different LOS, tripling the number of modes. Next, 
modes are grouped in spherical fc-shells with the width 
Ak/k = 0.186 (chosen as a trade-off between mini- 
mal mean x^ minimal numerical noise, but our 
results below are insensitive to this value). For each 
shell (the fc-dependence is implicit below), we shall 
fit the data set {/ti, P_^y(/ii), crp(/ii)} (where i runs 
through all modes within the shell) with the ansatz 
^atIa*) ~ X]j=i ^j-^jiP')^ where basis functions Xj{p) — 



{1,^ ,^, } and coefficients 
j — 1,2,3, respectively. 



{P^o,P^2,P^4}, for 



To minimize the merit func- 

2 



tion X = T.^[[PkTi^^i)-T.j=la]XJ{^l^)j /(jp{pi)\^ , 

we employ the standard General Linear Least Squares 
method (see, e.g. [HI). This results in best-fit coef- 
ficients Oj = ^jk Pk with associated error esti- 
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FIG. 1. Top left panel: the angle-averaged matter density power spectrum, Agg — Pssik) /2tt'^ (unitless), from N-body 

results for total matter (solid/black) and for /GM matter (green/dot), linear total density (dashed/blue), and 3PT total density 

- — -2 

(dot-dashed/magenta), at z = 6.757. Other panels: the total matter density power spectrum (multiplied by 5Ti,{z), the square 
of mean 21cm brightness temperature) as "measured" by fitting the mock data (i.e. reionization simulation) for the 3D power 
spectrum of 21cm brightness temperature anisotropy to the ^k-decomposition, A^4 = P^i{k) (in mK^), at different 
phases of reionization (mass- weighted global ionization fraction Xi^m = 0, 0.01, 0.10, 0.50, 0.68, 0.78, and 0.92, respectively). 
Red data points are reionization simulation results. Black data points assume homogeneously-ionized IGM with same global 
mean Xi^m as reionization simulation. Solid black curves are the expectation from linear scheme, i.e. equation Q evaluated 

- — -2 

using total density power spectra from N-body simulation (multiplied by STf,{z)). Fractional errors plotted in inset are with 
respect to the linear expectation, i.e. (p^is'-flt _ phncar-j^phnoar g^.j-Qj. bars and shaded regions correspond to the sampling 
variance of the simulation volume. 



mates <y{aj) — \fCjj- Here we define the 3x3 ma- 
trix ajk = Y^i Xj{\jii)Xu{pn) I a'\y{pii) , whose inverse is 
the covariance matrix C = and a 3- vector /3j — 

Mock data from a reionization simulation. We 

perform a new large-scale, high-resolution N-body sim- 
ulation of the ACDM universe (performed with the 
CubeP^M code [1, 0) in a comoving volume of ibox — 
425Mpc//i on each side using 5488^ (165 billion) par- 
ticles. To find halos, we use the spherical overdensity 
method and require them to consist of at least 20 N-body 
particles; this implies a minimum halo mass of 10^ Mq. 
We use subgrid modeling to compute the halo population 
with mass between 10* — 10^ M©. Assuming that the gas 
traces the CDM particles exactly, we grid the density and 
velocity in the IGM (i.e. halo mass excluded) on a 504^ 
grid using SPH-like smoothing with an adaptive kernel. 
Halo lists and density fields are then processed with the 
radiative transfer code C^Ray Each halo releases 
ionizing photons per baryon per At — 11.5 Myrs, with 
= 10 (/^ = 2) for halos below 10^ M© (above 10^ M©), 
respectively. To incorporate feedback from reionization, 
halos less massive than 10^ M© located in ionized regions 



do not produce any photons. We refer the readers to [9| 
for the elaboration of the codes and the feedback process, 
and Iliev et al. {in prep.) for the details of this simula- 
tion. The simulation used the following set of cosmolog- 
ical parameters I^a = 0.73, = 0.27, flh = 0.044, h = 
0.7, CTg = 0.8, ns = 0.96 where Hq = lOO/ikms"! Mpc^^, 
consistent with the seven-year results [loj . 

We then use the nonlinear density, velocity and ioniza- 
tion data from the simulation to calculate the nonlinear 
3D 21cm power spectrum, using the MM-RRM scheme 
in [3] for mapping N-body and reionization grid data in 
real space onto a redshift-space grid, and separate out 
the best-fit 4*'^-moment using the aforementioned angu- 
lar separation prescription. The linear and quasi-linear 
schemes both predict that this 4"'-moment should be 
given by equation ([3]), which we test by evaluating the 
r.h.s. of equation ^ directly from the simulation N-body 
results, for comparison. Some preliminary results were 
previously summarized by us in fll| . 

Results and discussions. In Figure [l] we plot the 
best-fit 4*''-moment of the fully nonlinear power spec- 
trum, and the benchmark linear expectation (i.e. equa- 
tion ^ evaluated with nonlinear simulation data). Note 
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FIG. 2. Systematic error of the separation scheme. We 
plot the fractional error of the best-fit 4*''-moment with re- 
spect to the linear expectation, in the inhomogeneous reion- 
ization case, as a function of global ionization fraction a;i,m, 
for k = 0.24,0.35, 0.52, 0.77, 0.94/1 Mpc"\ respectively. The 
solid/open dots correspond to the case where the systematic 
error is greater/less than the sampling error of our simulation. 
The inset is zoom-in to < a;; m < 0.10. 



that, while the mock 21cni signal is from the IGM, as 
is the observed signal, the total density power spectrum 
is expected in the linear scheme. (There is a difference 
between the total and IGM density power spectra, as in 
Figure [1] resulting from the exclusion of particle mass in 
halos when computing the IGM density.) We focus on the 
range of wavenumbers, 0.2 < k < 1/iMpc^^, for which 
k is large enough to avoid a large sampling variance but 
small enough to avoid a large aliasing effect. 

• Consistency check of the decomposition pipeline. We 
confirm that: (1) the total density power spectrum from 
N-body simulation agrees with the linear power spec- 
trum (from CAMB(T2) at large scales k < 0.5/iMpc~^ 
At smaller scales, it agrees with the result from third- 
order perturbation theory (3PT) , but is enhanced rel- 
ative to the linear power spectrum. (2) The best-fit 4***- 
moment agrees with the total density power spectrum 
at z = 30 when the IGM is everywhere neutral and the 
density fluctuations are of linear amplitude. 

• Effect of IGM nonlinear density /velocity fluctuations. 
For diagnostic purposes, we first investigate the case in 
which the ionized fraction at each point in the IGM is set 
equal to the (mass-weighted) global mean value Xi^m in 
the reionization simulation at that redshift. In this case, 
the best-fit 4*''-moment at different redshifts is enhanced 
with respect to the total density power spectrum, and the 
deviation increases from 20% [z ~ 14) to 70% [z ~ 7), 
as structure formation proceeds. These results show that 
nonlinear density and velocity fluctuations cause the /Zk- 
decomposition to make a systematic error, not caused by 
ionization patchiness. 

• Effect of inhomogeneous reionization and nonlinear ve- 



locity fluctuations. Early in reionization, the best-flt 
4*'^-moment for inhomogeneous reionization is suppressed 
relative to that for the homogeneously partially-ionized 
case. This is because fluctuations in neutral fraction 
and density anticorrelate in a universe reionized "inside- 
out," i.e. overdense regions are ionized earlier on aver- 
age than underdense regions. More precisely, inhomoge- 
neous reionization affects the anisotropy through its cou- 
pling with nonlinear velocity fluctuations, because reion- 
ization patchiness, alone, cannot introduce anisotropy 
in 21cm power spectrum (e.g. in the quasi-linear /Xk- 
decomposition scheme). A quantitative explanation will 
be formulated in more detail in Mao et al. {in prep.). 

This effect cancels the enhancement due to nonlinear 
fluctuations in IGM density and velocity, alone. Inci- 
dentally, the systematic error (i.e. departure of the best- 
fit 4"^-moment from linear expectation) first crosses zero 
(less than 10% at all scales) when Xi,„i « 10% (z ~ 9.6 in 
our simulation). Afterwards, this error grows to 60% for 
Xi^m ^5 50%. As determined by the competition between 
these two effects, this error depends on both the reion- 
ization epoch and the scale of interest — the smaller the 
ionized fraction Xi^m and the smaller the wavenumber fc, 
the smaller the error (see Figure [2]). 

As the typical size of ionized regions grows larger than 
the scales plotted here, the best-fit 4**^-moment at fc = 
0.5 — 1/iMpc^^ becomes less suppressed after the 50% 
ionized epoch. For this fc-range, the net error changes 
sign again when Xi^m — 68% (with error < 20%). 

At late epochs {xi^m ^ 0.8, z <7), the systematic er- 
ror for all scales is large, > 100%. This is due to the 
breakdown of the perturbative expansion, i.e. the expan- 
sion of the 3D 21cm power spectrum in neutral density 
fiuctuations becomes divergent when Sp^^^ > 0(1), as the 
ionized bubbles expand and percolate in the universe. In 
addition, the lightcone effect [l^ becomes non- negligible 
at this late time, and can further squeeze the anisotropic 
power spectrum along the LOS, as does the redshift-space 
distortion. Therefore, the estimate of 100% error here is 
only a lower limit to the actual error at late times. 

• The sampling variance. Our simulation volume is 
large enough that these systematic errors quoted above 
are all greater than the estimated sampling errors (except 
when the systematic error crosses zero), so they repre- 
sent statistically significant deviations from the expecta- 
tions of the /.tk-decomposition, rather than accidents of 
the particular numerical realization. 

Conclusion. This letter is the first attempt to quan- 
tify in detail the intrinsic precision of the linear scheme 
in which the cosmological matter density power spec- 
trum can be extracted from 21cm observations of the 
EOR. Two effects may spoil the extraction, a major one 
due to the coupling between inhomogeneous reionization 
and nonlinear velocity, and a minor one due to nonlinear 
density and velocity fluctuations alone. The competition 
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between these identifies two phases of reionization par- 
ticularly interesting to cosmology — Xi^m — 68%, where 
systematic error is within 20% for k = 0.5 — 1 hMpc~^, 
and Xi^m — 10%, where systematic error is within 10% 
for all wavenumbers. The epoch of exact crossover is 
likely to depend on the reionization scenario, as does the 
reionization patchiness. The 68%-ionized epoch is more 
accessible to first-generation observations (e.g. LOFAR, 
MWA 15]) because of smaller foreground and thermal 
noise at lower redshift. 

We summarize our results in Figure [D We see that, 
for ionized fraction less than 40% {z ~ 7.760), the lin- 
ear /ik-decomposition works well for large-scale mea- 
surement k < 0.24 /iMpc"^, with errors within 20%, 
and its application could still be acceptable at smaller 
scales, down to fc ~ 0.5/iMpc~^, with errors up to 
50%. During the intermediate phase of reionization 
(xi^m — 0.4 — 0.7), decomposition at the intermediate 
A:-range, k ~ 0.35 — 0.5 /i Mpc~^, is acceptable with er- 
rors up to 50%. However, in the late stage of reioniza- 
tion {xi^m ^ 0.8), it is difficult to use 21cm measurement 
to extract the cosmological information in the approach 
proposed by the linear scheme, because the systematic 
error is of order unity. Measurements at very large scale 
k < 0.2ft,Mpc~^ might still be useful at this late phase, 
but a survey volume much larger than that of LOFAR 
(which is comparable to our simulation volume) is neces- 
sary to reduce the sampling variance in the 4*^-moment 
significantly (e.g. MWA with 30° field of view). 
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